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DYNAMICAL EVOLUTION OF DISK GALAXIES 


By Frank Hohl 
Langley Research Center 

SUMMARY 

A computer model for isolated disks of stars is presented and is used to study the 
self-consistent motion of large numbers of point masses as they move in the plane of the 
galactic disk. The Control Data 6600 computer system at the Langley Research Center 
was used to integrate the equations of motion for each star for systems containing from 
50 000 to 200 000 stars. Any initially cold balanced disk was found to be violently 
unstable. A sufficient amount of velocity dispersion will stabilize all small-scale dis- 
turbances. However, most disks investigated were found to be unstable against slowly 
growing long -wavelength modes, and after about two rotations the disks tended to assume 
a bar-shaped structure. It was also found that the final mass distribution for most disks 
could be closely approximated by an exponential variation irrespective of the initial mass 
distribution. 

To study the development of spiral structure, the model was modified to include a 
fixed central force similar to that in the Schmidt model of the Galaxy. The mass of the 
stars in the disk was taken to be from 5 to 50 percent of the total mass of the Galaxy. 

The evolution of a number of initial distributions of stars was investigated. The results 
of the calculations gave a velocity dispersion for the disk stars which was about 50 per- 
cent larger than the value of about 30 km/sec found from observation of stars in the solar 
neighborhood. For some of the disks investigated, a pronounced spiral structure remained 
even after 8.5 rotations. 


INTRODUCTION 

A large fraction of the observable galaxies are spiral galaxies which have a disklike 
structure. De Vaucouleurs (ref. 1) finds that out of a sample of 1528 bright galaxies 
61.1 percent are normal spirals. Lin (ref. 2) estimates that about 70 percent of the known 
galaxies are normal spirals. A normal spiral galaxy consists of a relatively thin disk- 
shaped spiral pattern and a central nucleus. Although the disk can be treated as a two- 
dimensional system, the central nucleus is a complicated three -dimensional system. In 
general, it is found that the nucleus rotates like a solid body and that the disk rotates with 
differential rotation (ref. 3). Masses of galaxies range from 10^ to 10 12 solar masses 



and the dimensions are typically near 10 4 parsec (pc). The Galaxy has a mass of about 
10 solar masses. Most of the recent information on the interstellar matter and the 
structure of the Galaxy is contained in references 4 and 5. It appears that all spiral gal- 
axies contain a small amount of gas, probably less than 10 percent of the total mass. The 
mass of gas (mostly atomic hydrogen) in the Galaxy is about 3 percent of the total mass 
(ref. 6). There may be some additional gas in the form of molecular hydrogen (refs. 7 
and 8). Much of the gas is concentrated into a thin layer in the galactic disk between 
about 5 and 15 kpc (ref. 9). In the solar neighborhood, the gas constitutes about 15 
to 20 percent of the total mass in the spiral arms (refs. 10 and 11). In the spiral arms 
the gas density is about three times larger than in the interarm regions. Most of the 
bright young stars and H II (ionized hydrogen) regions are in the galactic plane near the 
spiral arms. It is found that many of the spiral galaxies show a two-arm spiral struc- 
ture. An excellent collection of photographs of galaxies has been compiled by Sandage 
(ref. 12). In that publication the detailed structure of a number of spiral galaxies is 
shown. 

Oort (ref. 13) has divided the problem concerning spiral structure into two parts. 
The first part is concerned with the origin of the spiral structure - namely, the mecha- 
nism that causes the spiral structure. The second part raises the question of the per- 
sistence of the spiral pattern. For example, the Galaxy is about 10 billion years old and 
its rotational period (in the solar region) is about 200 million years. Thus, a theory is 
needed which allows the spiral structure to persist even after 50 rotations of a galaxy. 

One difficulty in many of the early theories of spiral structure is that differential rotation 
will wind up any spiral structure in only a few rotations (ref. 14). 

Many of the recent efforts to explain spiral structure represent the spiral galaxy by 
an infinitesimally thin disk (refs. 15 to 19). In such models, the galaxy is considered to 
be an isolated system and the collective field due to purely gravitational effects deter- 
mines the motion of the stars. Lin and Shu (refs. 20 to 22) have used such a disk model 
to argue that the spiral structure persists because it is due to a density wave. Similar 
ideas were put forth earlier by Lindblad (refs. 23 to 25); however, his arguments were 
less convincing because they were supported primarily by his studies of individual star 
orbits. Recent investigations (refs. 26 and 27) show that the effects of finite thickness of 
the galactic disk do not change qualitatively the results obtained with the infinitesimally 
thin disk model. 

An assumption in the present report and in much of the recent work on spiral struc- 
ture is that the spiral pattern is caused primarily by gravitational effects. This assump- 
tion seems to be justified if a comparison is made of the energy density due to the differ- 
ent phenomena in the galaxy. For example, the ratio of the energy density in galactic 
rotation to that in the magnetic field (10-9 Wb/m2) is about 320 (ref. 22). The energy 
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density corresponding to cosmic rays, radiation, and turbulent gas motion is even smaller 
than that of the magnetic field. Nevertheless, the magnetic field may become important r 
(ref. 28) for the motion of gas clouds in the galactic disk (ref. 9). For example, 
Chandrasekhar and Fermi (ref. 29) proposed that the magnetic field may be responsible 
for the stability of the spiral structure. Greyber (ref. 30) assumes that matter falling 
inward from the galactic halo is guided into the nucleus by the magnetic dipole field of 
the galaxy. The gas is then ejected radially outward from the nucleus into the disk. 

Since analytical methods are difficult to apply to the dynamical evolution of galaxies, 
computer models have been developed to investigate the evolution of self -gravitating sys- 
tems. There are essentially two different types of computer models. In the first type of 
computer model, the binary interaction between any pair of stars is important during the 
evolution of the system under investigation and great care is taken to treat the near 
encounters of the stars very accurately. In fact, the formation of binaries is one of the 
difficulties in such calculations. Such a model applies to systems like globular clusters 
where collisional effects are important. The pioneering work in this area was done by 
Pasta and Ulam (ref. 31). Von Hoerner (refs. 32 and 33) made extensive calculations by 
following the three-dimensional motion of up to 25 stars. Aarseth (refs. 34 and 35) was 
able to treat the motion of 100 masses with his model for clusters of galaxies. More 
recently, R. Wielen (ref. 36) has used the model to perform a number of numerical exper- 
iments to study the dynamical evolution of star clusters. So far, the number of stars that 
have been treated with this model is only a few hundred because of the large amount of 
computer time required. 

The second type of computer model applies to systems that can be considered col- 
lisionless. The force acting on a particular star is then determined by the mean gravi- 
tational field of the system and the effect of nearby stars can be neglected (ref. 37, 
p. 363). Chandrasekhar (ref. 38) points out that for the Galaxy a star can describe at 
least 100 rotations about the galactic center before the effects of encounters with other 
nearby stars become important. A galaxy can therefore be studied by using a 
"collisionless" computer model in which close encounters are effectively neglected. The 
first such models used were simple one -dimensional models, where the stars are strati- 
fied into plane parallel layers or mass sheets (refs. 39 to 43). The one -dimensional 
model is useful as an approximation to the motion of stars normal to the galactic plane 
of the Galaxy. The one -dimensional model was also useful in investigating "violent 
relaxation" of self-gravitating systems (ref. 44). Henon (refs. 45 and 46) used a model 
of concentric spherical mass shells to study the dynamical mixing of spherical star 
clusters. A model of concentric mass rings was used by Toomre (ref. 18). Two- 
dimensional models, in which the motion of infinitely long mass rods is followed, have 
been used by Hockney (ref. 47) and by Hohl (refs. 48 to 50) in an attempt to investigate 
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galactic structure. The applicability of the two-dimensional rod model to stellar struc- 
ture is questionable. The rod-star model is likely to be valid only for very few galaxies, 
such as the cigar- shaped galaxy NGC 2685 (ref. 12). The results in the present report 
were obtained with a collisionless computer model which simulates the motion of large 
numbers of mass points or finite-length mass rods that are confined to move in the plane 
of the galactic disk. The model effectively simulates the evolution of an isolated disk of 
stars. Lindblad (refs. 51 and 52) pioneered such calculations by following the motion of 
up to 192 mutually attracting mass points in the given central field of the Galaxy. By 
placing the mass points initially in a system of concentric rings with circular velocities, 
Lindblad investigated the mutual disturbances in such a system to simulate the spiral 
structure of galaxies. Because Lindblad was able to follow the motion of only a rather 
small number of stars, the model is not effective for simulating the evolution of a colli- 
sionless galaxy. Miller and Prendergast (ref. 53) developed a model to study the motion 
of stars in a plane for systems which are doubly periodic and the forces, star positions, 
and velocities are allowed only discrete (integer) values which are less than some given 
maximum value. 


SYMBOLS 


D impact parameter 

E gravitational field, V<p 

G gravitational constant 

H Green' s function 

h displacement of star for epicyclic motion 

i = /3T 

i,j,k,l indices 

k wave number 

l rod length 

M total mass of disk, Nm 

m star mass 
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mass of field star 


mass of test star 

dimension of array of cells; total number of stars 

number of encounters 

star density 

potential energy 

pressure 

°r 
°r,min 

R 0 radius of disk 

r,0,z cylindrical coordinates 

T kinetic energy 

t time 

V rotational velocity 

velocity of test star before encounter 
Av£ increase in perpendicular velocity of test star 

x,y,z Cartesian coordinates 

T angular momentum 

6t time step for advancing motion of system 

6(z) Dirac delta function 


m£ 

m t 

N 

N c 

n 

P 

P 



k epicyclic frequency 

X wavelength 

X c largest unstable wavelength 



p mass density 

o surface mass density 

CT r ,Oq velocity dispersion in radial and azimuthal direction, respectively 

G r,min velocity dispersion defined by equation (45) 

t c relaxation time for collisional effects 

T r rotational period 

<p gravitational potential 

co angular velocity 

co Q angular velocity of cold balanced disk 

Subscripts: 

n,m cell location in x- and y-direction, respectively 

max maximum 

min minimum 

Notation: 

* scaled quantity 
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Fourier transformed quantity 


( ) average value 


COMPUTER MODEL 


The computer model used to investigate the dynamics of disk galaxies is illustrated 
in figure 1. The NxN array of cells shown in figure 1 is superposed over the plane of 
the galactic disk. The array of cells is introduced only for the purpose of calculating the 
gravitational potential. The cells are identified by n,m with n = 0,1,2, ...,N-1 
increasing in the x-direction and with m = 0,1,2,...,N-1 increasing in the y-direction. 
The cell in the lower left-hand corner is 0,0 and that in the upper right-hand corner of 
the array is N-1,N-1. The stars move over this imaginary array of cells. At the cen- 
ter of each cell a mass density is defined which is given by the number of stars in that 
cell. The number of stars in a cell is usually of the order of 100 and can become much 
larger near condensations. The density distribution is used to obtain the gravitational 
potential at the center of each cell. From the gravitational potential, the force acting at 
the position of a star is computed. Newton’s equations of motion are then used to advance 
the position and velocity of each star by a small time step. For the parameters of a 
typical galaxy, retardation or relativistic effects need not be considered. 

One complete cycle for advancing the motion of the system by a time fit consists 
of the following procedure. First, the distribution of mass a n m is used to obtain the 
gravitational potential cp n by effectively summing over the density. Second, the 
gravitational field at the position of the stars is computed from the potential cp n m . 

Third, by applying Newton's laws of motion, the motion of all the stars is advanced for a 
small time step fit. This procedure represents one cycle and it is repeated until the 
desired evolution of the system is achieved. 


Equations of Motion 

The motion of the stars is described by the differential equations 


and 


dVx _ 0_£ 

dt 9x 


dVy d (p 
dt 9y 


V x = 


dx 

dt 


= *L 
dt 


( 1 ) 


( 2 ) 
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The variable cp represents the gravitational potential and the gravitational field is 
given by E = V<p. For a star in the (n,m)th cell, equations (1) and (2) in the time- 
centered finite difference form (refs. 54 and 55) are 


H[ V x( t + f)- V x(‘-f)]=M n , m < 3) 

and 

£[x(t + 6t) - x(tj] = V x ^t + (4) 

with similar equations for the y -components. The numerical calculations can be speeded 
up greatly by scaling the distance so that the cell dimensions are equal to unity, that is, 
Ax* = Ay* =1. If, in addition, the velocity and mass of a star are scaled as 


= V. 


_5t_ 

Ax 


and 


2 (Ax) 2 


m 


( 5 ) 


( 6 ) 


then the potential is scaled as 


_ (5t) 
2 (Ax)2 


The equations of motion take on the simplified form 

< 7 > 


and 


x*(t + 6t) = x*(t) + V*(t + (8) 

Two methods were used to obtain the gravitational field. The simple method is to 
let each star in a particular cell experience the same field components, namely, 
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(®x)n,m ^Wl,m " ^n-l,m 


(9) 


for the x-component of the gravitational field in the n,m cell and 

( E y)n,m = ^m+l “ \m-l ^ 

for the y-component of the field. Equations (9) and (10) show that all stars in the cell 
(n,m) experience the same gravitational field and the value of the field will jump in 
crossing the cell boundaries. A smoother variation of the field acting on a star is 
obtained by means of a bilinear interpolation of the fields (as given by eqs. (9) and (10) at 
the four cell centers surrounding the position of a particular star. The two components 
of the field are then given by 

E* = (1 - 6 y)(l - 6x) (4 +l m - + «y(l - «x) (<4.i im+1 - <_!, m+1 ) 

+ 6x(l - 6y) (4. 2>m - < m ) + «x 6y - <_ m+1 ) (11) 

and 

E y - <1 - W - fc > «m+l - <111-1) + ^ 

+ 6y(l - «x) + 6x 5y (< +1>m+2 - (12) 

where the pertinent parameters are defined in figure 2. It is found that the bilinear 
interpolation gives a slightly more definite structure for the condensations which occurred 
during the initial evolution of a system. After about one rotation, the results obtained by 
the two methods display essentially the same structure. 

If a star should leave the NxN array of cells, the field acting on it is calculated 
by placing all the mass remaining in the system at the center of the array. The stars 
outside the array will not interact among themselves, but they will be attracted by the 
central force due to the mass placed at the center of the array. Whenever an appreciable 
number of stars leave the array, the calculations are no longer accurate and the computer 
run should be repeated by either increasing the array or by changing the initial conditions. 
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The Potential Calculation 


In calculations with the two-dimensional rod model (refs. 47 and 49), the gravita- 
tional potential is easily obtained by solving the two-dimensional Poisson equation 

+ Mr = 4irGp(x,y) (13) 

8x^ 8y^ 

An attempt therefore might be made to obtain the gravitational potential for the disk model 
by the same method. The Poisson equation then becomes 

2% + Mjr + = 47rGcr(x,y) 6(z) (14) 

8x^ dy* dz 6 

where 6(z) is the Dirac delta function and cr(x,y) is the surface density of stars in the 
plane of the disk. The difficulty is that no means are available to evaluate 8 %(p/dz%. 
Therefore, 


<p(x,y) = G C f g - -— - ■ ■ dx’ dy' (15) 

(/(x - x’) 2 + (y -y ’) 2 

is used to obtain the gravitational potential from the mass density (the primes denote 
variables over which integration is performed). Presently the density is given only at a 
finite number of cells separated by a unit distance so that the integral can be changed to 
a summation 




n,m 


N-l N-l 

■ 1 1 
i=0 j=0 


(16) 


where N is the dimension of the array of cells and H is Green’s function defined by 



1 


(17) 


To perform directly the summation indicated by equation (16) requires the summation of 
n 4 terms. For N = 100, = 10® and the time required to obtain cp becomes 

excessive. 
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Summation Method 


A faster approximate method to obtain the gravitational potential is illustrated in 
figure 3. In the summation method, the density cr n rn in each cell is first determined. 
Next, the density which consists of a n m plus the density of the surrounding 

eight cells is found. The potential at the center of a particular cell is then obtained by 
summing the contributions of the surrounding eight cells and the contribution of 0Vp at 
the center of each of the larger cells (fig. 3); that is, 

^m = pn,m + °h,m+l + °n,m-l + °h-l,m + °h+l,m + ^ (°n+l,m+l 

+ °h+l,m-l + °n-l,m-l + °h-l,m+l)j (^t)^ j**i-n,j-m 

i i 

where 

H i_ n i-m = ■ — - ((i - n) 2 + (j - m) 2 > o) (19a) 

ft i - n) 2 + (j - m)2 

H 00 = 0 (19b) 


The addition of the term <x_ w in equation (18) is equivalent to setting the self- 
potential equal to unity. The condition Hq q = 0 in equation (19b) simply omits the con- 
tribution of (<rJ in the double summation of equation (18). The values of Hi_ n i_ m 
V 1 /n,m ’ J 

are stored in an array and have to be calculated only once during a particular run. The 

indices i and j in. the double summation increase in increments of 3 such that one of 

the values of i and of j is equal to n and to m, respectively, and the summation is 

over every third value of (a T \ . With equation (18) the number of operations is 
_ _ ' ' n, m 

(N + 1) 2 [9 + (N + l) 2 / 9], a reduction of almost 1 order of magnitude over the number of 
operations required in the direct summation method. In a subsequent section the results 
obtained by the summation method are compared with those obtained by the Fourier 
method. In principle, the method could be extended to include a third class of cells con- 
taining the mass of 81 small cells o n m (or 9 large cells (Orp) V The number of 

9 r ’ V 9/ 1 K 

operations would then be (N + 1) [17 + (N + l)y 81J, a reduction of nearly 2 orders of 
magnitude over the number of operations in the direct summation method. 
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Fourier Method 


Another method to obtain the gravitational potential makes use of the fast Fourier 
transform methods now available (ref. 56). The Fourier transform of the density is 
defined as (ref. 57) 


°M = N 


N-l N-l 


n^O m=0 


£ °h,m ex P i |f( kn + lm ) 


Similarly, the Fourier transform of Green's function is given by 

„ N-l N-l 


S k,l=(|) I £ H n>m exp[?I(lm + lm)] 
' n=0 m=0 J 


Applying the finite convolution theorem to equation (16) gives the result (refs. 57, 58, 
and 59) 


N-l N-l 


Z Z °i,j H i-n,j-i 


i=0 j=0 


N-l N-l 


* (I) £ X 5 k,l g k,l “Pp 1 If 1 + lm) ] 

' ' k=0 1=0 


From equation (22) and the definition of the Fourier transform, it is clear that 


^k,l = \lK,l 

Therefore, the potential (p n>m is obtained directly from the inverse Fourier transform 
of a • H. Such a method gives a doubly periodic system and one can expect the tidal 
interaction between image galaxies to be important. A similar method is used by Miller 
and Prendergast (ref. 53) in their investigation of doubly periodic stellar systems. 

The Fourier transform method just described can be modified to obtain the potential 
distribution for an isolated system. This modification is achieved by increasing the num- 
ber of cells by a factor of 4 and by confining the system to one-quarter of the array of 
cells. The mass density in the remaining three-quarters of the array will then always 
be identically zero. 

Consider now that in addition to the array under investigation, the summation in 
equation (16) is extended over all the doubly infinite array of images. However, Green's 
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function H n?m is now modified so that it corresponds to the correct single particle 
potential for particle separation r less than N/2 (one -half the dimension of the array) 
and to zero interaction for r greater than N/2. Even though the system is still doubly 
periodic, there is no longer any interaction between adjacent image systems because 
their masses are separated by at least N/2. 

Thus, to get the correct potential for an isolated system at the expense of a four- 
fold increase in storage, Green’s function to be used in equation (22) is 

(°Sn,m S f) 

H n ,m = f ' (23a) 

yn^ + m^ (n^ + m^ / OJ 


®N-n,m **n,N-m ®N-n,N-m ®n,m 


(23b) 


and 


H 


0,0 


= 1 


(23c) 


H 0,0 =1 

The Fourier transform of 


As before, setting 
unity. 

symmetry of H 


is equivalent to setting the self -potential of a star equal to 
H n m need be done only once. Also, because of the 

'N 


n,m 


only a finite cosine transform on a ^ + lj x + lj 


mesh is 


required. The modified Fourier transform approach is described in references 58 
and 59. It should be pointed out that the Fourier transform method solves equation (16) 
for the isolated system exactly (within computer rounding error). A computer listing and 
a description of the Fourier method are given in the appendix. 


Arbitrary Force Law 

The two methods presented for obtaining the gravitational potential can easily be 
extended to three-dimensional problems. Also, the force law between particles can 
easily be changed by simply changing Green's function H. For example, some of the 
effects of finite thickness of the galactic disk can be simulated by using finite -length mass 
rods instead of point masses to represent the stars. The force of attraction F between 
two mass rods of unit mass and of length l alined perpendicular to the galactic plane 
with their centers in the galactic plane is 

r 


F = 


2G 

i 2 


vM^)- 
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where r is the separation of the two rods. Green's function then becomes 



where r nm = y/n^ + Figure 4 shows H nin for four values of l. 

Comparison of Methods for Obtaining Gravitation Potential 

The time required to obtain the gravitational potential by the Fourier method for 
a 64 x 64 active array (Fourier transform must be performed on 128 x 128 array) was 
2 seconds on the Control Data 6600 computer system at the Langley Research Center. 

The Fourier analysis and synthesis were performed with a program written in COMPASS 
assembly code. The remainder of the potential solver was written in FORTRAN IV. For 
the summation method, the time required was 6 seconds with the program written in 
COMPASS assembly code. The time required for a FORTRAN IV program is increased 
by nearly a factor of 3 . 

The time required to advance the motion and to calculate the density for a 
50 000-star system is 8.6 seconds when the positions and velocities of the stars are 
stored on magnetic disks and are placed into fast storage 10 000 at a time. 

Since the star coordinates are stored on magnetic disks, the number of stars that 
can be handled is essentially limited only by the available computer time. With packing, 
the summation method only takes one -half of the storage required by the Fourier method. 
Since the potential array is always kept in fast storage, the summation method makes 
possible the use of a larger mesh. 


Energy and Momentum Conservation 
The angular momentum of the disk at time t is given by 



where mj is the mass of the ith star, V x ^ is the x-component of the velocity of that 
star, and the summation extends over all the stars. For most of the systems investi- 
gated, the angular momentum remained constant to within 0.15 percent during the first 
one -half rotation of the system. 
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The kinetic energy of the disk at time t is 



An approximate expression for the potential energy is 


P ~ "15/ 2 °h,m^n,m 
n m 


( 26 ) 


(27) 


A better definition of the potential energy should be devised by using the definition of 
potential energy as the work done on a test star and apply it to the present model. The 
difficulty in the use of equation (27) is that in regions of large mass condensations this 
equation gives a value of <p which is too large for the potential energy. During the first 
one-half rotation, equations (26) and (27) give a value of the total energy which remains 
constant within 0.5 percent. 


Effect of Varying Discretization Parameters 

Figure 5 shows how the evolution of the system is affected by varying the discreti- 
zation parameters of the system. This figure illustrates the evolution of initially uni- 
formly rotating (solid-body rotation) balanced disks of stars. In figure 5(b) the initially 
cold disk contains 50 000 stars and the gravitational potential for this disk is obtained 
by the summation method. The time indicated is in rotational periods. Figure 5(a) pre- 
sents the evolution of an identical system but with the potential calculated by means of 
the Fourier transform technique. It should be pointed out that the initial perturbations 
present in the systems in figures 5(a) and 5(b) are identical since the initial positions for 
the two disks were generated by the same pseudorandom number generator. The number 
of time steps per rotation is 200. After a quarter rotation, the summation method is 
seen to give a finer structure for the condensations than the Fourier transform method. 
The distance between condensations in figure 5(b) is about the dimension of one large cell 
(3x3 small cells) and the condensations are within one large cell. Because of the 
approximation made by introducing the larger cells in the summation method, the inter- 
action between masses separated by two (or more) cell dimensions is decreased. This 
effect explains the finer structure of the initial condensation in figure 5(b). A comparison 
of figure 5(a) and figure 5(b) shows that after one-half rotation the overall structure of 
the condensations is the same for the two methods. Even after one rotation the results 
are very similar. It should be noted that since the system is violently unstable, any 
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small change in method of computation or initial conditions will cause rather large devia- 
tion in the overall structure of the condensations which occur as the disk breaks up. 

Figure 5(c) shows the effect of changing the time step from 200 to 400 time steps 
per rotation with the gravitational potential being obtained by the Fourier method. For 
the first one-half rotation the evolution is the same as that of figure 5(a). After one rota- 
tion the number of condensations for the two disks is the same, but the relative positions 
are changed somewhat. Again, it should be noted that for all disks in figure 5, the initial 
pseudorandom positions are identical so that the perturbations caused by them are the 
same. 

Figure 5(d) shows the effect of increasing the number of stars to 200 000. This 
increase is obtained by distributing four stars in the neighborhood of the position of each 
of the original 50 000 stars. Thus, the perturbations caused by the initial positions of the 
stars should be similar to those of the previous systems. The results shown in figure 5(d) 
indicate that increasing the number of stars has very little effect on the evolution of the 
system. The effect of reducing the mesh size was also checked. It was found that on a 
128 x 128 active mesh the condensations are much finer than those for the 64 x 64 mesh. 
This result is to be expected for the extreme case of an initially cold disk. Toomre 
(ref. 18) has shown that a cold disk of stars is unstable to perturbations of arbitrarily 
small wavelengths. Since the smallest condensations possible in the present model are 
of the dimensions of a cell, a decrease in the cell dimensions will allow smaller conden- 
sations to occur. It has been shown by Toomre (ref. 18) that the short -wavelength modes 
are stabilized by the effects of velocity dispersion. The effect of changing the cell dimen- 
sion was also compared for an initial disk of stars with sufficiently large random veloc- 
ities to stabilize modes of the order of a few cell dimensions. Decreasing the cell dimen- 
sion by a factor of 2 had essentially no effect on the evolution of the system (ref. 59). 

It was also found that using the cloud-in-cell (CIC) (ref. 60) method had little effect 
on the evolution of the system. In the CIC method of calculation, each star has the dimen- 
sion of a cell and its mass is distributed among four cells in proportion to the fraction of 
particle area in a cell. 

The photographs in figure 5 of the star disks were made with the model dd80B 
display-recorder system of the Data Display Division of Control Data Corporation. This 
system produces both 3 5 -mm motion -picture film or page -size recordings. 

EVOLUTION OF SELF-CONSISTENT DISK GALAXIES 

The first disk galaxies to be investigated are completely self-consistent disks of 
stars without any imposed gravitational field. The initial conditions chosen are such that 
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the gravitational attraction acting on the stars is balanced by the centrifugal force due to 
the rotation of the stars around the center of the system. 


Initial Condition 

Various authors (refs. 17 and 61 to 63) have constructed time -independent models 
for disk galaxies. An initial condition used extensively in the present calculations corre- 
sponds to the familiar simple model of a uniformly rotating disk with a surface mass 
density given by 

a(r) = cr(0)^l -^J /2 




(28) 


where r is the radial coordinate, R 0 is the radius of the disk, and M is the total 
mass of the disk. The uniform angular velocity required to balance the gravitational 
attraction is then 


co 


= JsGirtim/lIt* 


(29) 


where M = Nm, N being the number of stars in the system and m the mass per star. 
To treat disks with a higher central condensation, Toomre suggested the use of the fol- 
lowing generalization of equation (28): 

.. / o \ n-1 / 2 

a n (r) = -M^(2n+ 1) 1 - M (30) 

2»S| \ R 2 j 

where n = 1, 2, 3, . . , . 

Dr. C. Hunter of the Massachusetts Institute of Technology has found that the 
angular velocity co for such a density distribution is given by 


where 


co 


2 [ t \ _ 7tGM 


©9 = 


2 4n -v.) 2 E"-i)0 v 



(31) 


(32) 
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and 2 F i> the usual hypergeometric function (ref. 64), is found from 



The square of the rotation speed is 


Vg(r) = r 2 o> 2 (£) (34) 

For the first five values of n, equation (34) has the form 



(35a) 

v 2 <r,= ^(^) 2 (i)( 1+3{2 ) 

(35b) 

v i< r) = T(i) (25e)( 1 + 2{2 + 554 ) 

(35c) 


(35d) 

V 5 (r) ” + T« 2 + T« 4 + 4 « 6 + 9 « 8 ) 

(35e) 


Figure 6 shows cr n (r) as a function of r/R Q for the first five values of n. The corre- 
sponding curves for the values of the rotational velocity are given in figure 7. The initial 
density distribution for most of the disks of stars investigated corresponds to n = 1 in 
equation (30); a few disks with n = 3 and n = 5 were also studied. The initial positions 
of the stars are obtained by using a pseudorandom number generator which gives a uni- 
form distribution between 0 and 1 to generate three numbers x r , y r , and z r . If 
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(36) 


xf + y? < 1 


and 

(l - x 2 - Yp) > z r (37) 

the initial position of a star is given by the x,y coordinates 

x(t = 0) = 2R 0 (x r - 0.5) (38) 

and 

y(t = 0) = 2R 0 (y r - 0.5) (39) 

The initial circular velocity of the stars is given by equation (35a). If the stars have an 
initial velocity dispersion, a uniform velocity dispersion cr r = is superposed on the 
circular velocity given by equation (35a). For some of the uniformly rotating disks, the 
pressure effects due to the initial velocity dispersion were taken into account in deter- 
mining the initial angular velocity required to balance the disk. For these disks, the 
velocity dispersion was taken to be a function of radius so that 

) f^ m 

The angular velocity required to initially balance the warm disk for <r r = a q is given by 


w 


2 — w 2 + 
o 


1 8p 
rcr 8r 


= co 2 
o 




R; 


(41) 


The angular velocity u> 0 is that of the cold disk given by equation (29); the pressure p 
is found from 

p = a(r)a 2 

where the density a is given by equation (28). If <r r (r) ^ cr 0 (r), the angular velocity 
required to balance the disk is 
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(42) 


co 2 = co 2 + 


ra(r) 9r 


|a(r)a 2 (r)] + -^-jo^(r) - o 2 ^)] 


The initial circular velocity required to balance a given density distribution, which 
is a function of r only, can be obtained directly from the computer model. The given 
initial density distribution is used to calculate the gravitational potential and the radial 
gravitational field. The radial field E r (r) is then used to obtain the angular velocity 
for stars at a radius r by equating the centrifugal acceleration to the gravitational 
attraction 


or 




= E 


’( r )| 


(43) 



(44) 


since 


V 0 (r) = rco(r) 

Effect of Velocity Dispersion 

All the calculations for uniformly rotating disks presented in this section were per- 
formed with disks containing 50 000 stars. There are 200 time steps per rotation and 
the gravitational field was calculated by equations (9) and (10). A 64 x 64 active mesh 
was used for the potential calculation. 

By using a dimensional order-of-magnitude argument, Toomre (ref. 18) predicted 
that a cold (zero velocity dispersion) infinitesimally thin disk of stars is violently 
unstable. The growth time of the unstable modes is approximately proportional to the 
square root of their wavelength (ref. 16) and, therefore,, the small -wavelength distur- 
bances will be the most unstable. Toomre (ref. 18) finds that the root-mean-square 
radial velocity dispersion required to stabilize all axisymmetric disturbances anywhere 
in the disk is 


a r,min 


3 


.36 — 

K 


(45) 


where a and k are the local values of the density and epicyclic frequency, respec- 
tively. The epicyclic frequency is given by 
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( 46 ) 


k 2 = i_£ + 3 9£ 
a r 2 r 9r 

Initially the cold uniformly rotating disk is balanced and 

2 

8£_ 3irGNm r 

9r r 4 t?2 


2 2 

where Vg(r) = V-^(r) is given by equation (35a). Substituting equation (47) into equa- 
tion (46) gives the following equation for k %: 


3irGNm 2V 0 W 


The expressions for the density cr(r) and (eqs. (28) and (48), respectively) can be 
used in equation (45) to obtain the radial velocity dispersion needed to stabilize the disk. 
Thus, 

a r,mln = Wl ‘ [0 ‘ 49 > 

The longest axisymmetric unstable wavelength of a thin pressure free disk was 
given by Toomre as 

\ _ 47r 2 Gcr / c ^ 


For the uniformly rotating disk, equation (50) reduces to 


~ Rot / 1 



Figure 8 shows the evolution of an initially uniformly rotating balanced disk with 
zero velocity dispersion. The time is in rotational periods given by 


T r = 277 


37iGNm 
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The initial kinetic energy of the cold balanced disk is 


v 0 

T(t=0) = J ( cr(r)Va(r)2vT dr 
2 Jq 0 

_ 3ttGM 2 
20R o 


2 2 

where a(r) and V|(r) = V|(r) are given by equations (28) and (35a), respectively. The 
initial potential energy is given by 


P(t=0) 


1 r R ° 

2 A) 


cr(r)(p(r)2irr 


dr 


37rGM 2 

ior 0 


where 


a>M = 37rGM r 2 - 37rGM 
8R 3 4R o 

The virial theorem is of course satisfied at t = 0, that is, 

2T(t=0) + P(t=0) = 0 

The values of P(t=0) and T(t=0) were found to agree with the computer generated 
values obtained by using equations (26) and (27). As predicted by Toomre (ref. 18) the 
cold disk is violently unstable and the growth time for the instabilities is a fraction of 
the rotational period. Since the growth time of a disturbance is approximately propor- 
tional to the square root of its wavelength, the fast appearance of the small-scale con- 
densations is to be expected. After 1.1 rotations there remain five main condensations. 
The same information is contained in figure 9 which shows the evolution of the gravita- 
tional potential contours in the plane of the disk. The relative sizes of the condensations 
can be more easily distinguished in figure 9. The evolution displayed in figures 8 and 9 
is rather violent and the stars quickly build up a large velocity dispersion. The random- 
ization can also be observed in the evolution of the velocity distribution for this system 
(fig. 10). At t = 0.4r r , the velocity distribution displays a ringlike structure which, 
however, has disappeared after 1.1 rotations. The five smaller clusters which remained 
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after 1.1 rotations in figure 8 exist mainly because the pressure caused by the increased 
velocity dispersion of the stars counteracts the gravitational attraction towards the center 
of each cluster. If the trajectories of individual stars are plotted, it is found that they 
oscillate in the potential wells set up by the five condensations of stars. This fast ran- 
domization observed in figures 8 and 10 must be caused by collective effects or violent 
relaxation as discussed by Lynden-Bell (refs. 44 and 65). 


The relaxation time for collisional effects r c can be estimated by summing the 
effects of the encounters of individual stars. This summation can be performed by 
applying the test star approach (ref. 38) to a plane stellar system. Consider a test star 
of mass m^ which encounters a field star of mass mj . Let the impact parameter and 
the initial relative velocity be D and V, respectively. The orbits of the two stars are 
hyperbolas with asymptotes that make an angle it - 20 in the center-of-gravity coordi- 
nate system, where 


sin 0 = 



D 2 V 4 


1 - 1/2 


G 2 (m t + m f y 


(53) 


The angular deflection of the test star is therefore 20. Before the encounter, the test 
star has a velocity 


v, = 2 — V 

t m^ + m^ 

in the center-of-mass coordinate system. The increase in the perpendicular velocity of 
the test star caused by the encounter is given by 


( 4v i) 


2 _ t, 2 erfr.2 


vf sin' 5 20 


4m 


- — - V 2 sin 2 0 cos 2 0 


From equation (53) it is found that 

PV/o 2 (mf+m t f 


COS ' 5 0 


1 + 


D 2 V 4 


G 2 (m f + m t y 


(54) 
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Equation (54) then becomes 



4m|v 6 D 2 


C 2 (m f + m t y 


1 + 


G 2 ( 


D 2 vf 

m f + m t ) 5 


So far, the derivation is the familiar one found in many textbooks. The next step is to 
sum over encounters in the plane of the disk galaxy rather than sum over a volume. The 
number of encounters N C (D) per unit time with impact parameter between D and 
D + dD is 


N C (D) dD = 2nV dD 

where n is the number density of field stars. The total increase in the transverse 
velocity is then given by 

.min 

SGVmfn / V^in/G^ + m t ) V^rnax/^ + m t) 


nu + m* 


i2 < 1 + D 


mm 


-|2' 


G (m f + m t ) 


2/1 + D 


max 
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G^mj + m^ 


+ i-tan-1 
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4G 2 m 2 n 
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V^Dmax 

- I tan" 1 

V 2 D • 
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or 
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and 


D max >:> D min 


where D m in is the cell dimension and D max corresponds to the dimension of the sys- 
tem. The time required for a 90° deflection corresponds to vj* = or 


T c = 


V 3 Pmin 

l2m? r 


^v|^> 4G‘ s m|n 


(55) 


The angular velocity of the uniformly rotating disk is given by equation (29) and 
initially the mean quadratic difference of velocity between two stars in the disk is 

-srf 


(56) 


Also, the rotational period of the uniformly rotating disk is 

v l/2 

2ir 


4R; 


Tr « * n \.3irGNm/ 


■ = 2tt 


sV 


(57) 


The ratio of the collisional relaxation time r c to the rotational period r r is then 



2N * 2 x 10"“ N 


(58) 


where n is expressed as N/irR^ and D m i n ^R 0 * The estimate for the collisional 

relaxation time t c given by equation (58) shows that for a 50 000-star disk, collisional 
effects are important only after about 1000 rotations. 

The velocity given by equation (56) is probably too large, and a better estimate for 
the velocity to be used in equation (55) is the minimum velocity dispersion required for 
stabilizing the disk as calculated by Toomre (ref. 18) 


V = 3.36 — 


The relaxation time then becomes 


where m = ^ is used. 


r c * 9.5 
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mm 


/c^m 


(59) 
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With the following approximations 


and 


2w 2 « /c 2 


7tRqGct 2 


ffR o ff 

m 


« N 


equation (59) can be written 


0.32 



( 60 ) 


The ratio of the relaxation time to the rotational period (for r r = 2 t t/w) is 

— = 0.052 * 2 X 10" 3 N (61) 

and for N = 50 000 the system can be considered collisionless for about 100 rotations. 

Toomre (ref. 18) predicted that the short-wavelength modes are stabilized by the 
effects of velocity dispersion. To check his predictions, all stars in the initially uni- 
formly rotating disk were given a uniform velocity dispersion. Figure 11 shows the 
evolution for such a disk of stars where the root-mean-square value of the random veloc- 
ity is 0.068V^(R o ^ or 6.8 percent of the circular velocity at the edge of the cold balanced 
disk. The initial positions of the stars are the same as those for the disk of stars in fig- 
ure 8. The disk is still violently unstable but, when compared with the evolution of the 
cold disk in figure 8, it can be seen that the initial condensations take longer to form and 
are no longer as concentrated as those for the cold disk. The smallest condensations 
formed are now of the order of three cell dimensions. 

The evolution of the disk with a velocity dispersion equal to 13.6 percent of the 
circular velocity at the edge of the cold disk is shown in figure 12. The smallest con- 
densations formed are now of the order of six cell dimensions. The evolution of the disk 
shows that the small -wavelength modes have been stabilized. However, the disk is still 
violently unstable with respect to large -wavelength modes. 
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An initial velocity dispersion of 20.4 percent of the circular velocity at the edge of 
the cold balanced disk results in the evolution displayed in figure 13. The disk is now 
nearly stabilized, but it still forms a barred spirallike condensation about 10 cell dimen- 
sions across. Figure 14 shows the corresponding evolution in velocity space wherein the 
system quickly attains a ringlike distribution. 

Finally, the velocity dispersion was increased to 27.2 percent of the circular veloc- 
ity at the edge of the cold balanced disk. The resulting evolution presented in figure 15 
shows that the overall system is now relatively stable. Equation (49) gives the minimum 
velocity dispersion to stabilize the disk at r = 0 as 34.2 percent of the circular velocity 
of the cold balanced disk. The average velocity dispersion required to stabilize the disk 
should be less than the maximum value. Therefore, the value of 27.2 percent found from 
the computer simulation is in good agreement with the predictions of Toomre (ref. 18). 
Figure 15 indicates that an increase in the central density occurs as the system evolves. 
This effect is illustrated in figure 16 where the average density is plotted as a function of 
radius during the evolution of the system. The density was obtained from the number of 
stars with radii between r and r + Ar. 

For the evolution of the disks in figures 11 to 15, the velocity dispersion was super- 
posed on the circular velocity of the stars for the cold balanced disk. Therefore, the 
disks were not exactly balanced at t = 0, since the added velocity dispersion results in 
a pressure which causes the disk to expand. Also, the velocity dispersion was isotropic - 
that is, a r = = Constant. By using a velocity dispersion a r (r) as given by equa- 

tion (40) and a velocity dispersion Og = o v = o r (for the present disk) the system 

£U?q 

has the velocity dispersion required by Toomre (ref. 18) for stability. The disk is ini- 
tially balanced by reducing the initial solid-body rotation according to equation (41). Fig- 
ure 17 presents the evolution of such a balanced disk for which the initial velocity disper- 
sion is one -half that given by equation (49). As can be seen, the system is still unstable 
and quickly forms a number of condensations. The evolution is similar to that in fig- 
ure 12. For a velocity dispersion equal to that given by equation (49), the evolution of 
the initially balanced disk is as shown in figure 18. The only structure which slowly 
develops is a mode with dimensions of the size of the whole disk. These long -wavelength 
disturbances have an azimuthal dependence and are not covered by the theory of Toomre. 

For the results given in figures 8 to 18, the gravitational potential was obtained by 
using the Fourier transform method; also, the evolution is shown only up to 1.1 rotations. 
In figure 19 is presented the evolution of an initially balanced cold disk of 50 000 stars 
when the potential is obtained by the summation method. The evolution is shown up to 
2.2 rotations. The long time evolution (from 1.4 to 2.2 rotations) is typical for unstable 
disks. For disks with slightly different initial conditions, it has been found that after two 
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to five rotations the system consisted of one large central condensation around which 
orbited one to three small condensations. Initially uniformly rotating disks that are 
stabilized against small-scale disturbances (such as the system in fig. 18) tend to assume 
a barlike structure after about two rotations. 

Effect of Finite Thickness 

Some of the effects of finite thickness of the disk of stars can be taken into consid- 
eration by representing the stars as finite-length mass rods. The main effect in using 
finite -length rod stars is to decrease the close range interaction for the stars as shown 
in figure 4. The long-range collective effects are essentially unchanged. The initial 
density distribution of mass rods used to investigate the effects of finite thickness is 
taken to be that given by equation (28). The total mass of each rod is m. Since the 
gravitational field resulting from a distribution of mass rods is different from that due 
to a distribution of point masses, the rotation required to balance the cold disk is no 
longer given by equation (29). Therefore, the gravitational field is calculated by the 
model from the given initial density and equation (44) is used to obtain the initial rotation. 
Figure 20 shows the evolution of a 50 000-mass-rod system; each rod has a length of two 
cell dimensions. The structures displayed in figures 8 and 20 are very similar to 
t = 0.3r r to t = 0.5r r . After t = 0.5r r , the slightly different initial conditions cause 
rather large deviations in the overall structure of the condensations which occur as the 
disks break up. These deviations are to be expected since the systems are violently 
unstable. The evolution of the gravitational potential in the plane of the disk is presented 
in figure 21. Note that the gravitational potential is scaled differently from that in fig- 
ure 9. Figure 22 shows the circular velocities of the stars as a function of radius for 
the rod system in figure 20. At t = 0, it can be seen that the circular velocity increases 
nearly linearly with radius, except near the edge of the disk where the circular velocity 
becomes constant. This leveling off of the circular velocity of course will not occur 
when the initial velocity is obtained by using equation (29). The evolution given in fig- 
ure 22 indicates that the stars quickly increase their random velocities. Figure 23 shows 
how the distribution of the radial velocities evolves in time. Initially all the stars have 
zero radial velocities. As the system breaks up, the stars gain radial velocities and 
move along elliptical trajectories in the V r ,r space. 

When the length of the mass rods is increased to five cell dimensions, the effect of 
finite thickness becomes more noticeable as is shown in figure 24. The short-range 
interaction is now reduced appreciably and, as a result, the growth time for small-scale 
instabilities is decreased. Also, as can be seen from figure 24, the number of conden- 
sations which form initially has been reduced when compared with those for the zero- 
thickness disk (fig. 8). Figure 25 shows the evolution of the radial velocities of the 
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stars for the system in figure 24. Because of the decrease in the short-range inter- 
actions of the mass rods, the radial (or random) velocities increase at a much slower 
rate than for the point mass system. It appears that the main effect of increasing the 
length of the mass rods is to slow down the growth rate of the small-scale disturbances. 

Effect of Increased Central Condensation 

Hunter (ref. 15) finds that increasing the central condensation of a self -gravitating 
disk decreases the growth rate of some of the unstable axisymmetric oscillations. To 
investigate the effect of an increased central density, an initially cold balanced disk of 
stars with a density a^(r), given by equation (30), is studied. The evolution of this 
stellar system is shown in figure 26. As expected, the evolution displayed in the figure 
indicates that an increased central condensation alone does not stabilize the system 
against large-scale disturbances. The time is given in units normalized to the rotation 
periods of the cold uniformly rotating disk. As can be seen from figure 7, the period 
corresponding to stars at the edge of the disk is then about 25 percent larger than that 
corresponding to the normalized time given in figure 26. 

The effect of superposing a uniform velocity dispersion was also investigated and 
the results are shown in figures 27 and 28. In figure 27, the uniform velocity dispersion 
is 30 percent of the circular velocity at the edge of the cold balanced disk (see fig. 7 
for n = 3). Most of the small-scale condensations no longer form, and after one rota- 
tion the system displays an S-shape structure. If the uniform velocity dispersion is 
increased to 60 percent of the circular velocity at the edge of the cold balanced disk, the 
evolution displayed in figure 28 is obtained. Apparently the system initially cools by 
ejecting the hot stars and it then contracts. After t = 0.7 r r , the system expands again. 

The evolution of a system with an even higher central condensation given by Og(r) 
(fig. 6) is presented in figure 29. The system is initially cold and balanced and the time 
shown is normalized to the period of a uniformly rotating balanced disk with the same 
total mass and radius. This system is violently unstable and quickly condenses into a 
circular ring as suggested by Hunter (ref. 16). The ring-shaped condensation later 
breaks up into four smaller condensations. Figure 30 shows the evolution of the same 
system, but with the stars replaced by mass rods each of which has a length of five cell 
dimensions. The system no longer condenses into a circular ring. After t = 2.2r r , 
most of the mass is concentrated in the large central condensation. 

Gaussian Mass Distribution 

One of the problems of great interest is that of a galaxy being perturbed by the 
close passage of some other galaxy. An example of particular interest is the proposed 
close passage of the large Magellanic Cloud at 20 kpc from the center of the Galaxy some 
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500 million years ago (ref. 66). However, before such a problem can be undertaken by 
the present methods, a stable disk galaxy must first be generated which can then be per- 
turbed by the passage of the Magellanic Cloud. As a first attempt at a stable system, a 
disk with a Gaussian mass distribution is considered. The Gaussian mass distribution 
is given by 

°g( r ) = cr(0)exp £-7! r r 2j (62) 

where M is the total mass of the disk 

M = Og(r)27rrdr 

and a(0) is the mass density at r = 0. For the present calculations, the value of 
7rcr(0) /M was taken to be 2 x 10"^. The rotational velocity required to balance a cold 
disk with a Gaussian mass distribution has been calculated by Toomre (ref. 17) as 
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is the confluent hypergeometric function with a,p,z as arbitrary variables. For an 
initially cold balanced disk with a Gaussian mass distribution, the evolution is similar 
to that shown in figure 29; this was to be expected since for large n,<7 n approaches a 
Gaussian distribution. To determine whether Toomre' s minimum velocity dispersion 
(ref. 18) stabilizes a disk with a Gaussian mass distribution, the initial velocity disper- 
sion was taken to be 
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where a r min and x(r) are given by equations (45) and (46), respectively, and 
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Because of the added velocity dispersion, the initial angular velocity of the stars co(r) 
for a balanced disk is lower than co 0 (r) and is given by equation (42). Figure 31 shows 
the variation of u> 0 , c o, and k with r for a disk of stars with a Gaussian mass dis- 
tribution. The evolution of the disk is given in figure 32. There are 50 000 stars in the 
disk and 200 time steps per rotation. The time is given in units of the rotational period 
of the cold balanced disk 2n/<x> 0 (r) at r = 10 kpc. The rectangular border enclosing 
the disk is at x = ±19 kpc and y = ±19 kpc. The evolution displayed in figure 32 indi- 
cates that even though the disk is stabilized against axisymmetric disturbances according 
to Toomre’s formula, the system is still violently unstable against nonaxisymmetric dis- 
turbances. After 1.5 rotations, the system has assumed a bar-shaped structure which 
varies very little during the following two rotations. The change (error) in the total 
energy of the system is less than 1 percent for the three rotations and the change in the 
angular momentum is about 0.1 percent. It should also be noted that at t = 0 the virial 
theorem is satisfied, that is, 2T(t=0) = -P(t=0). Figure 33 shows the variation of the 
radial velocities of the stars for the stellar system in figure 32. It can be seen that the 
velocity dispersion of the stars increases with time and that the stars tend to concentrate 
in regions of small r (as is obvious from fig. 32). The corresponding evolution in 
V x ,Vy space is shown in figure 34. In order to get more quantitative information than 
can be obtained from figures 32 to 34, the disk was divided into a number of concentric 
rings each with a width of 1/2 kpc. The radial dependence of various parameters aver- 
aged azimuthally over each ring was then obtained. Figure 35 presents the radial velocity 
dispersion obtained in this manner for the disk of stars in figure 32. Initially, the velocity 
dispersion for larger r increases very rapidly as can be seen from a comparison of the 
results at t=0 and t=0.25r r . Again, the time shown is in units of 2ir/w 0 at 
r = 10 kpc. After three rotations, the radial velocity dispersion in the central part of the 
disk has increased to about 170 km/sec from the initial value of 105 km/ sec. The vari- 
ation of Q = a r /a r>m i n is given in figure 36. The value of cr r min is obtained from 
the azimuthally averaged parameters of the disk at time t by finite difference methods. 
After the system has settled down at t = 3T r , Q has a value of approximately 4 over 
most of the disk. Because of this large value of Q, no spiral structure is expected to 
appear. 

In interpreting the results shown in figures 35 and 36 the nonaxisymmetric shape 
of the disk at certain times should be kept in mind. Figure 37 presents the azimuthally 
averaged mass density as a function of radius. These results show that the central mass 
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density of the disk increased by a factor of 4 during the first two rotations. The dash- 
line curve shown at t = 3r r corresponds to an exponential density distribution with a 
scale length of 1.5 kpc and was obtained from 

cr(r) cc exp(-r/1.5) 

It can be seen that the exponential variation closely fits the mass distribution of the disk. 
The variation of the radial gravitational field for the disk with an initial Gaussian distri- 
bution is shown in figure 38. The gravitational field is obtained by averaging the field 
for x between 0 and -15 kpc and x between 0 and +15 kpc. It should be noted that the 
gravitational field at t = 3r r is much closer to that of the Galaxy (ref. 4) than is the 
gravitational field for the initial Gaussian mass distribution. Four individual star orbits 
taken at random from the 50 000 stars in the disk are shown in figure 39. The orbits 
indicate that stars initially near the center of the disk have a tendency to become trapped 
in even tighter orbits as the central mass density increases. Stars farther out have a 
tendency to escape from the system. To give an indication of the time step used in the 
integration, orbit 4 in figure 39 is plotted as a sequence of dots, one dot at the end of 
each time step. 

The foregoing results indicate that the disk of stars with a Gaussian mass distribu- 
tion and with a given initial velocity dispersion as required to stabilize all axisymmetric 
modes (ref. 18) is not stable. Increasing the initial velocity dispersion by 10 percent and 
using a different random number sequence for the initial star positions did not change the 
evolution of the disk very much. The result can be seen in figure 40 which shows the 
evolution of such a disk. The effect of using finite -length mass rods instead of point 
masses is illustrated in figure 41 with an initially balanced disk having a Gaussian mass 
distribution and zero velocity dispersion. The length of the mass rods is 2 kpc (4 cell 
dimensions). The time shown is in units of 27r/u> 0 (r) at r = 10 kpc for the thin disk, 
that is, the same units used in figures 32 and 40. It is seen in figure 41 that the small- 
scale modes are stabilized by the reduction of the interaction potential for stars in close 
proximity. However, large-scale modes are still violently unstable and at t = 2.5T r 
the system has assumed a bar -shaped structure. Note that the rectangular border 
enclosing each frame in figure 41 is at x = ±17 kpc and y = ±17 kpc. The results 
shown in figures 32 to 41 indicate that a disk of stars with a Gaussian mass distribution 
is not a stable configuration. 


Exponential Mass Distribution 

For the disk of stars with an initial Gaussian mass distribution, it was found that 
the final mass density was closely described by an exponential distribution. It is 
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therefore of interest to investigate whether a disk of stars with an initial exponential 
mass distribution will remain stable. For the disk investigated, the initial mass density 
is given by 


a(r) = cr(0)exp (-r/r Q ) 

where cr(0) is the central density and r 0 = 3 kpc. Figure 42 shows the initial values 
of k , u> 0 , and w as a function of r as obtained by finite difference methods from the 
computer model. The evolution of the disk is shown in figure 43. The initial mass den- 
sity was cut off at 15 kpc (or five scale lengths). This cutoff should not affect the 
dynamics of the disk since only 4 percent of the total mass for the exponential distribu- 
tion lies outside 15 kpc. The time step used in the calculations is very small (800 time 
steps per rotation at a radius of 10 kpc). The evolution of this disk is not quite as violent 
as the evolution of the disk with the Gaussian mass distribution. (Compare figs. 43 and 
and 32.) However, after only 1.5 rotations the system has also assumed a bar-shaped 
mass distribution with a halo of stars moving in larger orbits around the central bar. 

Only slow changes in the shape of the system occur after t = 1.5r r . The variation of the 
velocities is shown in figure 44. A more quantitative view of the radial velocity disper- 
sion can be obtained from figure 45 where cr r is plotted as a function of r at different 
times during the evolution. At t = 1.5T r , the radial velocity dispersion has increased 
by only about 50 percent over the initial value. The variation of Q = o - r /^ r ,min * s 
shown in figure 46. Note that Q increases slowly to a value of about 2; compare this 
value with that of about 4 obtained for the disk with a Gaussian mass density (fig. 36). 
From the mass distribution shown in figure 47, a tendency for the central mass concen- 
tration to increase can be observed. In fact, the final mass density is closely described 
by an exponential distribution given by 

a(r) °c exp(-r/1.5) 

as indicated by the dash-line curve at t = 1.5r r . This is the same final density varia- 
tion as was found for the disk with a Gaussian mass distribution. Many of the disk 
galaxies investigated approach an exponential mass distribution; this fact is significant, 
especially since there is observational evidence (ref. 3) that the luminosity distribution 
of many spiral galaxies can be approximated by an exponential distribution. The varia- 
tion of the radial gravitational field with radius is shown in figure 48. As is to be 
expected from the increase in the central mass density, the field increases sharply near 
the center of the disk. The epicylic frequency k can be obtained directly from the 
following equation: 
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K 2 = 
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where E r = dcp/dr. The variation of the epicyclic frequency with radius is presented 
in figure 49. As expected, the epicyclic frequency increases appreciably in the central 
part of the disk. The orbits of six of the 50 000 stars in the disk are shown in figure 50. 
Note that the initial rotation is in the clockwise direction. As was observed for the 
Gaussian mass distribution, stars with small initial radii tend to be trapped near the 
center; whereas stars with large initial radii tend to move farther away from the center. 

It is of interest to compare the evolution of the "stabilized" exponential disk in fig- 
ure 43 with the evolution of an initially cold (zero velocity dispersion) disk. Figure 51 
shows the evolution of a cold disk of stars with the same initial exponential mass distri- 
bution as the disk in figure 43. The disk now shows many fine scale condensations and 
it quickly breaks up into a number of separate clusters of stars. Further work is being 
done on the exponential disk with the hope that a stable axisymmetric disk can be obtained. 


Comparison of Disk Model With Infinite -Mass -Rod Model 

It is of interest to compare the results obtained for a disk of point masses with the 
results obtained with cylindrical Systems that consist of infinitely long mass rods 
(refs. 47 to 50). Figure 52 shows the evolution of an initially nonrotating cold disk of 
stars. The time is in units of the rotation period of a balanced uniformly rotating disk. 
The system is seen to quickly collapse, then expand again, and present some filamentary 
structure. After t = 1.0r r , most of the mass has condensed into one large cluster of 
stars at the center of the system. These results are essentially the same as those 
obtained for the cylindrical galaxy (fig. 5 of ref. 49). 

The evolution of a cold disk with an initial rotation equal to one -half that required 
to balance the disk is shown in figure 53. Again, the evolution is very similar to that 
obtained for infinitely long mass rods (fig. 7 of ref. 49). The main difference is that for 
the point masses any structure is formed much quicker than for the rod stars. This 
stronger interaction is to be expected from the two-star interaction potential shown in 
figure 4 as a function of rod length. 

For the cylindrical galaxy (infinitely long rods), it was found that the uniformly 
rotating cylinder was unstable (ref. 50). By letting half of the stars rotate clockwise and 
half rotate counterclockwise, a stable system was obtained (fig. 16 of ref. 50). Figure 54 
shows the corresponding evolution of an initially balanced disk of point stars where 
25 000 stars rotate clockwise and the remaining 25 000 stars rotate counterclockwise. 
The system is violently unstable and quickly condenses into a ring shape. Thereafter, 
large random velocities quickly build up and cause the system to thermalize. 
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EVOLUTION OF SPIRAL STRUCTURE 


To study the development of spiral structure, the computer model was modified to 
include a fixed central potential in addition to the self-consistent disk population of stars. 
The central field is taken to represent the halo population and the central core of the 
Galaxy, which appear to produce a relatively time -independent field. The two central 
potentials used correspond to circular velocities (km/sec) given by 


. _ 4500r 
6 81 + r 2 


(64) 


and 
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y + =■ (65) 
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In figure 55 the two velocity curves are compared with the Schmidt model (ref. 4) of the 
Galaxy. The evolution of a system is affected little by using one or the other of the two 
rotation curves. The mass corresponding to the fixed central potential was taken to be 
as small as 50 percent of the total mass of the system. However, for most of the calcu- 
lations the self-consistent disk population contained from 10 percent to 20 percent of the 
mass of the system. 

Figure 56 gives the evolution of a 50 000-star system with a fixed central potential 
corresponding to the rotation curve given by equation (64). The disk stars represent 
10 percent of the total mass of the system. The time step used for this calculation is 
rather large, there being only 80 time steps per rotation at 10 kpc from the center of the 
system. Initially, the disk stars have a density of the form given by equation (28) with a 
disk radius of 20 kpc, zero radial velocities, and rotational velocities just sufficient to 
balance the gravitational attraction. It can be seen from figure 56 that the system quickly 
develops a multiarm spiral structure. The time shown in figure 56 and in all subsequent 
figures is in units of the rotational period at a radius of 10 kpc. The spiral structure 
remains even after 8.5 rotations, at which time the calculation was terminated. From a 
purely kinematical viewpoint, differential rotation should have wound up the spiral struc- 
ture before 8.5 rotations. In many of the calculations, there was a tendency to form cir- 
cular rings of stars in the central region of the disk. Figure 56 shows that after 2.5 rota- 
tions the system develops a number of concentric rings of stars which appear to be quite 
stable. In the outer parts of the system, the rings open up into a spiral structure. After 
6.5 rotations the system resembles a two-arm spiral galaxy with a tightly wound spiral 
structure in the central part of the system. The circular velocities of the stars for this 
galactic system are plotted as a function of distance from the center of the system in 
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figure 57. The velocity dispersion builds up substantially in the region above 10 kpc but 
remains quite small near the center of the system. This effect can be explained by the 
rather low value of the self-consistent gravitational field in the central region of the 
galaxy. 

To investigate more closely the appearance and stability of the concentric rings 
which appeared in figure 56, the same galactic system was investigated again but with an 
added initial perturbation. Equations (7) and (8) indicate that to start the calculation cor- 
rectly with the given velocity of a star at time zero V(t=0), the position of the star must 
be given at time t = 0 + 6t/2 or r(t = 6t/2). The perturbation was introduced by initially 
placing the stars at position r(t=0) with velocity V(t=0); this caused a shift of the cor- 
rect initial position of each star given by 


where 


Ar = r(t=0) 


{■ 



w(r) 


V fl (r) 


r 


( 66 ) 


and 


r = x + iy 


The perturbation is essentially a rotation (around the disk center) of the initial star posi- 
tions, about 2° from their equilibrium positions. The resulting motion of the stars in the 
plane of the disk is shown in figure 58. The rings are now much more pronounced and 
are quite stable. Again, the rings open up into a spiral structure in the outer parts of 
the disk. The rings are actually density waves which start at the center of the system 
and move radially outward. No net mass is carried radially outward by these density 
waves. The circular velocities of the stars for this galactic system are plotted as a 
function of radius in figure 59. The wavelike pattern of the distribution is caused by the 
epicyclic motion of the stars in the central field. This effect can be seen even more 
clearly in figure 60 wherein the radial velocities of the stars are plotted as a function of 
radius. The number of maxima and minima corresponds to the concentric rings in fig- 
ure 58. This behavior can be explained in terms of the epicyclic motion of the individual 
stars from the following equation describing the displacement of a star h(r,t) from its 
nearly circular orbit (ref. 38): 
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h(r,t) = r cos k( r)t 


(67) 



For the present system, the epicyclic frequency /c( r) is approximately 
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where equation (46) was used with 
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and Vq is given by equation (64). The radial wave number is 


k(r,t) = 



54 x 4500r 

(81 + r2) 5/2 


t 


and increases linearly with time. This is exactly the behavior indicated in figures 58 
to 60 where the distance between concentric rings or the wavelength given by 


X(r,t) 


2t r _ tt(81 + r 2 ) 5/2 
k(r,t) 27 x 4500rt 


( 68 ) 


decreases with increasing time. 

Observational evidence indicates that a ring or system of rings is a feature of 
many galaxies. Examples of spiral galaxies which display a structure of concentric 
circular rings of stars, similar to that shown in figure 58, are NGC 488 and NGC 1398 
(ref. 12). 

Another 50 000- star system with a fixed central potential given by equation (65) was 
investigated and the evolution of the system is shown in figure 61. The disk stars con- 
tain 20 percent of the total mass of the galaxy and the initial radius of the disk is 15 kpc. 
The time step used is such that initially there are 200 time steps per rotation at 10 kpc. 
The system quickly develops spiral structure which remains quite pronounced up to about 
four rotations. After about six rotations, the spiral structure becomes quite diffuse due 
to the buildup of random velocities. The buildup of the random velocities can be seen in 
figure 62, where the radial velocities of the stars are plotted as a function of the star 
radius. Initially, the disk is cold and all stars have zero radial velocities. As the evolu- 
tion of the system continues, the random velocities of the stars continue to increase up to 
about two rotations; thereafter, there is practically no increase. The same effect can be 
seen in figure 63, where the circular velocities of the stars are plotted as a function of 
the star radius. Figures 62 and 63 are plotted to the same scale. It should be noted that 
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in the central part of the system the self-consistent gravitational field due to the disk 
stars is rather small compared with the fixed central field. This fact explains the low 
value of the velocity dispersion in the central region of the disk (figs. 62 and 63). To 
study some of the individual star orbits more closely, the orbits of four stars have been 
plotted in figure 64. The motion of the stars is in the clockwise direction. To determine 
whether the stars followed an epicyclic motion, the radius of the star at time t was 
subtracted from the initial radius and the difference was plotted as a function of time in 
figure 65. The oscillations shown occur at the epicyclic frequency and the amplitude of 
the oscillations increases with increasing radius of the star orbit. The variation of the 
mass density with radius is shown in figure 66. The density is obtained from the number 
of stars with radius between r and r + Ar. In the absence of the applied central field, 
the central increase in density would be more pronounced. Figure 67 shows the varia- 
tion of the average root-mean- square radial velocities with radius. Again, after two 
rotations there is little further increase in the velocity dispersion. The velocity disper- 
sion in the azimuthal direction is found to be smaller than the radial velocity dispersion 
by a factor given by 
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(69) 


This result is in agreement with theory. The minimum velocity dispersion required to 
stabilize all axisymmetric disturbances (ref. 18) 


<7 r> min=3.36S2 


is plotted in figure 68 for various times during the evolution of the system. The values 
of a and k are calculated from the parameters of the system at the given times. 

The variation of 


°r,min 

with r is plotted in figure 69. The value of Q appears to settle down to about 1.7. 
This rather large value of Q would, of course, suppress any further instabilities which 
could result in spiral structure. The diffuse, rather structureless appearance of the 
system shown in figure 61 at t = 6.8T r is to be expected because of the large value 
of Q. Toomre (ref. 18) calculates the longest unstable wavelength in the absense of 
velocity dispersion as 
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A plot of X c as a function of radius for various times during the evolution of the system 
is presented in figure 70. It can be seen that in the central part of the system the longest 
unstable wavelength becomes very small. This result is in agreement with the results 
shown in figure 61, where only very small-scale structure is visible in the central part 
of the system. 

The evolution of a disk galaxy with a fixed central field corresponding to equa- 
tion (64) and with 50 000 disk stars containing 10 percent of the total mass of the system 
is shown in figure 71. The initial density distribution is given by equation (28). The 
disk is initially cold and balanced with a radius of 15 kpc and there are 100 time steps 
per rotation at 10 kpc. The evolution of the system is similar to that shown in figure 61. 
At first, a multiarm spiral structure appears in the outer regions of the disk. As the 
evolution proceeds, the increase in the velocity dispersion causes any structure to 
become quite diffuse. Four star orbits from the stellar system in figure 71 and the 
deviations of these orbits from circular orbits are given in figure 72 and figure 73, 
respectively. The oscillations shown in figure 73 occur at the epicyclic frequency. The 
variation of the velocity dispersion with radius for this stellar system is presented in 
figure 74. For r = 0 to about 5 kpc, the dynamics is dominated by the fixed central 

field and the velocity dispersion remains quite low in this region. After three rotations 

oy» 

there is little further increase in the velocity dispersion. The variation Q = 

°r,min 

with r for this system is shown in figure 75. After three rotations, Q is well above 1 
for radii larger than about 5 kpc. 

The four systems investigated (figs. 56, 58, 61, and 71) had very little initial central 
condensation of the disk stars. To determine the effect of increasing the central mass 
density, the evolution of a disk of stars with an initial Gaussian mass density was studied 
(fig. 76). The 50 000 disk stars contain 20 percent of the total mass of the galaxy. There 
is now much more fine scale structure in the central portion of the disk when compared 
with structure in the other four systems. However, this structure is difficult to detect in 
figure 76. Abetter indication of the fine scale structure can be obtained from the evolu- 
tion of the same system in V x ,Vy space (fig. 77). The individual orbits of four stars 
in the system are shown in figure 78. The corresponding deviations of the star motion 
from circular orbits is shown in figure 79. The oscillations occur at the epicyclic fre- 
quency which remains practically unchanged during the evolution as can be seen from 
figure 80. The variation of the velocity dispersion with r is given in figure 81. The 
increased central star density has the effect of increasing the velocity dispersion in the 
central region of the disk. After about two rotations the velocity dispersion increases 
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very little in the outer regions of the disk. Because in the central part of the disk the 
fixed field is still large compared with the self-consistent field, the velocity dispersion 
in that region takes longer to reach a larger value. The variation of Q with r for 
this disk of stars is shown in figure 82. From the low value of Q for r ^ 5 kpc and 
from the increased central density (compared with the previous 4 systems), one could 
expect the fine scale structure which persisted in the central part of the disk. The den- 
sity in the center of the disk did increase by about 30 percent during the evolution shown 
in figures 56 to 62. 

Also investigated was a disk of stars with the same initial conditions as those for 
the system in figure 76, except that the mass of the 50 000 disk stars was increased to 
50 percent of the total mass of the system. The evolution of that system is shown in 
figure 83. The disk is initially violently unstable but quickly approaches a rather steady 
state with a much increased central density. By t = 1.75r r , the density at the center of 
the disk has increased by a factor of 6. The variation of the velocity dispersion with r 
is shown in figure 84. The large central density of disk stars now causes a rapid 
increase in the velocity dispersion in the central region of the disk. At t = 1.75r r , the 
value of Q is about 2 for the region 0 g r § 13 and it increases to about 4 for larger 
radii. The value of Q equal to about 4 indicates that after less than two rotations the 
system has essentially reached a rather hot steady state. The variation of the radial 
gravitational field with r is shown in figure 85. As the central mass density of the disk 
stars increases, the gravitational field due to the disk stars increases sharply near the 
center of the disk. 


CONCLUDING REMARKS 

A computer model for isolated disk galaxies has been described in detail. The 
model was tested by varying the discretization parameters. The evolution of the system 
was affected very little for the range of parameters used. The evolution of a cold bal- 
anced disk was calculated and the disk was found to be violently unstable; this result was 
in agreement with theoretical predictions. Adding increasing amounts of velocity dis- 
persion to the uniformly rotating disk had the effect of increasing the smallest size of the 
star condensations that were formed. A velocity dispersion equal to about 27 percent of 
the circular velocity at the edge of the cold balanced disk stabilized the disk against 
breakup and any local condensations. These results are in general agreement with the 
theoretical predictions of Toomre. However, the disks were found to be unstable against 
long -wavelength nonaxisymmetric disturbances and after only about two rotations the 
disks assumed a bar-shaped structure. Similar results are obtained for disks with a 
Gaussian density distribution and with a velocity dispersion equal to or even greater than 
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the "stabilizing" velocity dispersion calculated by Toomre. A disk with an initial expo- 
nential density distribution and with Toomre' s stabilizing velocity dispersion displayed 
the least violent evolution even though it also took on a bar shape. One interesting result 
of the calculations is that the final density distribution for the disk galaxies investigated 
can be closely approximated by an exponential distribution irrespective of the initial 
conditions. This result is in agreement with observational evidence which indicates that 
the luminosity distribution of many spiral galaxies can be approximated by an exponential 
distribution. To study the evolution of spiral structure, the model was modified by adding 
a fixed central force similar to that of the Schmidt model of the galaxy. The fixed force 
represents the halo population of stars and the central core of the Galaxy. For most of 
the systems with the fixed central force which were investigated, a spiral structure 
quickly appeared but became quite <liffuse after about six rotations. Also, the velocity 
dispersion of the disk stars increased to a value that was approximately 50 to 100 percent 
larger than the value of about 30 km/ sec for stars in the solar neighborhood. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., April 13, 1970. 
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APPENDIX 


COMPUTER PROGRAM FOR OBTAINING GRAVITATIONAL 
POTENTIAL OF ISOLATED DISK GALAXIES 


The following listing is a FORTRAN version of the subroutine GETPHI which is 
used to obtain the gravitational potential for isolated disk galaxies. 


SUBROUTINE GETPH I 

COMMON Z( 1025) ,Y< 1025) . RHO < 1 28 « 1 28 > • I2A, I TEST * S 

D I MENS I ON G ( 65 * 65 ) 

IF ( I TEST »EQ» 0 ) GO TO 40 
I2B=I2A-1 

N=2*+I2A 

N02=N/2 

M=N+1 

N21=N02+1 

NP2*N+2 

NGSQ=N21*N21 

IC=N02/2+l 

RNI»l . / ( N*N ) 

BS = f*J=l «N2 1 
DO 1 1 = 1 «N2l 

IF ( I *EQ.l .AND* J.EQ. 1) GO TO 1 
R2= U-l • )*( 1-1 • )+( J-I • )*( J-l • ) 

R=SQRT { R2 ) 

SS=SQRT< 1 .+S2/R2) 

G ( I » J ) = RNI* (2./S)*( IR/S)* ( 1 .-SS l+ALOG I.S/R+SS > ) 

1 CONT I NUE~ 

G ( 1 » 1 ) = G(1 *2 ) 

CALL GETSETC2, I2B) 

DO 3 J = 1 «N21 
DO 4 1*1 «N21 
4 Z( I >*G( I « J) 

CALL FTRANS (2. I2B) 

DO 10 1 = 1 « N2 1 
10 G ( I « J ) * Y ( 1 ) 

3 CONTINUE 

DO 2 1=1 » N2 1 
DO 19 J=1 ,N21 
19 Z( J)=G( I « J) 

CALL FTRANS (2» I2B) 

DO 22 J= 1 « N2 1 
22 G ( I » J ) *Y ( J ) 

2 CONTINUE 
40 CONTINUE 

CALL 6ETSET ( 3 « I 2A ) 

DO 7 J=1 «N 
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DO 8 1=1 *N 

8 Z< I ) *RHO ( I * J ) 

CALL FTRANS ( 3. I2A) 

DO 9 I - 1 « N 

9 RHO ( I « J)=Y< I) 

7 CONTINUE 

DO 25 1=1 «N 
DO 26 J=1 *N 

26 Z ( J ) --RHO ( I * J ) 

CALL FTRANSO, I2A ) 

DO 27 J=1 «N 

27 RHO ( I * J)=Y< J ) 

25 CONTINUE 

DO 11 J=2,N02 
DO 11 I =2 « NG2 
I 1=I+N02 
J1 =J+N02 

RHO ( I * J ) =RHO (I ♦ J ) *G ( I * J ) 

RHO ( 1 1 ♦ J ) =RHO ( 1 1 » J)*G( I *J) 

RHO ( I « J1 ) * RHO ( I * J1 )*G< I ♦ J) 

RHO ( 1 1 • J1 ) =RHO ( 1 1 « J 1 ) *G ( I ♦ J ) 

11 CONTINUE 

DO 37 I =2 * N02 

RHO ( I +N02 * 1 ) =RHO ( I +N02 « 1 ) *G ( I « 1 ) 

RHO ( I+N02»N21 >=RHO< I+N02»N21 > *G ( I «N21 ) 
RHO (1*1) =RHO ( I « 1 > *G ( I * 1 ) 

37 RHO ( I *N21 )=RHO( I»N21 )*G ( I «N21 ) 

DO 38 J=2 «N02 

RHOd * J ) =RHO ( 1 » J )*G ( 1 * J ) 

RHO ( 1 « J+N02 )=RHO ( 1 . J+N02 >*G ( 1 * J) 

RHO ( N2l »J)=RH0(N21*J)*G(N21 »J) 

38 RHO (N21 * J+N02)=RH0(N21 « J+N02 ) *6 ( N21 * J) 
RHOd « 1 )»RHO (1 * 1 >*G< 1*1) 

RHO ( N21 , 1 ) =RHO (N21 , 1 )*G (N21 « 1 ) 

RHOd * N21 )=RHO( 1 «N21 )*G ( 1 *N21 ) 

RHO (N21 *N21 )=RH0(N21 *N21 )*G(N21 *N21 ) 
CALL GETSET(4, I2A ) 

DO 14 J=1 *N 
DO 12 1=1 «N 

12 Z ( I ) =RHO ( 1 * J ) 

CALL FTRANS (4*1 2 A ) 

DO 13 I =1 *N 

13 RHO ( I * J ) =Y< 1 ) 

14 CONTINUE 

DO 28 I = 1 ♦ N 
DO 29 J=1 *N 

29 Z ( J ) =RHO ( I *J) 

CALL FTRANS (4* I2A) 

DO 30 J=1 *N 

30 RHO ( I » J ) =Y ( J ) 

28 CONTINUE 

I TEST = 0 
RETURN 
END 
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The interaction potential G(I,J) calculated in the listing is for finite -length mass 
rods as given by equation (24). The common variable S defines the length of the mass 
rods. The variable I2A defines the size of the rectangular array RHO(I,J) used in the 
potential calculations. When the subroutine GETPHI is called, RHO(I,J) contains the 
mass density and GETPHI places the values of the corresponding gravitational potential 
in RHO(I,J). The subroutine FTRANS(I,I2B) has been written by R. Hockney (ref. 58) 
and it performs a finite Fourier analysis or synthesis on the common input array Z 
and puts the result in the common output array Y. The subroutine performs a cosine 
analysis or synthesis for I = 2, a periodic analysis for 1=3, and a periodic synthesis 
for 1=4. The subroutine GETSET(I,I2B) initializes FTRANS and is called every time 
the arguments of FTRANS (I, I2B) are changed. 
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Figure 1.- Computer model illustrating the N x N array of cells used in calculating the gravitational potential. 
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(a) Cold disk containing 50000 stars,- 200 time steps per rotation; potential calculated by Fourier method on a 64x 64 active array of cells, 

(b) Cold disk containing 50000 stars; 200 time steps per rotation; potential calculated by summation method on a 64x64 active array of cells 

(c) Cold disk containing 50000 stars; 400 time steps per rotation; potential calculated by Fourier method on a 64x 64 active array of cells, 

(d) Cold disk containing 200000 stars,- 200 time steps per rotation,- potential calculated by Fourier method 

on a 64x 64 active array of cells. 

Figure 5.- Effect of varying the model discretization parameters on the evolution of a cold disk of stars. (Time in rotational periods given 

by eq. (52'.) 
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Figure 12.- Evolution of a disk of stars with a constant initial velocity dispersion equal to 13.6 percent of the circular velocity at the edge of 

the cold balanced disk. (Time in rotational periods of the cold disk.) 
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Figure 13.- Evolution of a disk of stars with a constant initial velocity dispersion equal to 20.4 percent of the circular velocity at the edge of 

the cold balanced disk. (Time in rotational periods of the cold disk.) 
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Figure 15.- Evolution of a disk of stars with a constant initial velocity dispersion equal to 27.2 percent of the circular velocity at the edge 

of the cold balanced disk. (Time in rotational periods of the cold disk.) 
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Figure 16.- Azimuthally averaged radial density for the system in figure 15. (Time in rotational periods of the cold disk.) 
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t = 0.3 t = 0.4 t = 0.5 t = 0.3 t = 0.4 t = 0.5 




t = 0.9 t = 1.0 t = l.l t = 0.9 t = 1.0 t = l.l 


Figure 27.- Evolution of a disk of stars with initial conditions the same as those 
of the disk in figure 20 except that a uniform velocity dispersion equal to 
30 percent of the velocity at the edge of the cold balanced disk has been added. 
(Time in rotational periods of the cold disk.) 


Figure 26.- Evolution of an initially balanced cold disk of stars with a density 
distribution given by o 3 (r). (Time in rotational periods of the cold disk.) 









t = 0.9 t = 1.0 t = l.l t = 0.9 t = 1.0 t = 1.1 


Figure 28.- Evolution of a disk of stars with initial conditions the same as Figure 29.- Evolution of an initially cold balanced disk of stars with a density 

those of the disk In figure 26 except that a uniform velocity dispersion equal distribution given by a 5 (r). (Time in rotational periods of the cold disk.) 

to 60 percent of the velocity at the edge of the cold balanced disk has been 
added. (Time in rotational periods of the cold disk.) 


















Figure 30.- Evolution of the stellar system displayed in figure 29 except that 
the point stars are mass rods having a length of five cell dimensions. 
(Time in rotational periods of the cold disk.) 


Figure 31.- Variation of o) Qt to, and k with r for a disk of stars with a 
Gaussian mass density. 
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Figure 36.- Variation of Q with radius for the Gaussian disk. {Time in rotational periods of the cold balanced Gaussian disk at r = 10 kpc.) 
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Figure 41.- Evolution of an initially balanced cold Gaussian disk consisting of Figure 42.- Variation of <y <o t and k with radius for a disk of stars with an 

mass rods each 2 kpc long. (Time in rotational periods of the cold balanced exponential mass distribution. 

Gaussian disk at r = 10 kpc.) 












Figure 45.- Variation of radial velocity dispersion with radius for the exponential disk. (Time in rotational periods of the cold exponential 

disk at r = 10 kpc.) 
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Figure 46.- Variation of Q = o r /o r m j n with radius for the exponential disk. (Time in rotational periods of the cold exponential disk 

at r = 10 kpc.) 
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Figure 47.- Variation of mass density with radius for the exponential disk. (Time in rotational periods of the cold exponential disk 

at r = 10 kpc.) 
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Figure 48.- 
























t = 1.0 


i.x 


i = 


al rotation equal to one-half that required to balance the disk. (Time in 
periods given by eq. (52).) 




















oas/ui: 


A Schmidt model 

V * 4600r/(100 + r 2 ) + 210r/(0.25 + r 2 ) 


■V fl = 4500r/(81 + r ) 



Figure 55.- Comparison of two rotation curves with the Schmidt model of the Galaxy. A fixed gravitational potential corresponding to the 

two rotation curves is used to study the evolution of spiral structure. 
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Figure 56.- Evolution of a 50 000-star system under influence of a fixed central force corresponding to rotation curve given by equation (64). 
The disk stars contain 10 percent of total mass of system. (Time in rotational periods of the initial balanced disk at r = 10 kpc.) 
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t = 2.500 t = 5.000 t = 8.500 

Figure 57.- Evolution of the circular velocities of the stars for the galaxy in figure 56. (Time in rotational periods of the initial balanced 

disk at r = 10 kpc.) 


90 












Figure 58.- Evolution of the same stellar system displayed in figure 56 except that the initial star positions have received an initial perturbation 
which is proportional to the distance of the star from the center of the disk. (Time in rotational periods of the initial balanced disk at 
r = 10 kpc.) 
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Figure 61.- Evolution of a disk of stars under the influence of a fixed central force corresponding to the rotation curve given by equation (65). 
The 50 000 disk stars contain 20 percent of the total mass of the system. (Time in rotational periods of the initial balanced disk at 
r = 10 kpc.) 
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Figure 71.- Evolution of an initially balanced cold .disk of stars under the influence of a fixed central force corresponding to the rotation curve 
given by equation (64). The 50 000 disk stars contain 10 percent of total mass of the galaxy. (Time in rotational periods of the initial balanced 
disk at r = 10 kpc.) 
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Figure 76.- Evolution of a disk of stars with a Gaussian mass distribution under the influence of a fixed radial force corresponding to the 
rotation curve given by equation (64). The 50000 disk stars contain 20 percent of the total mass of the galaxy. (Time in rotational periods 
of the initial balanced disk at r = 10 kpc.) 
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Figure 82.- Variation of Q = a r /a r m j n with radius for the system in figure 76. (Time in rotational periods of the initial balanced disk 

at r = 10 kpc.) 







Figure 83.- Evolution of the stellar system displayed in figure 76 except that the 500(H) disk stars contain 50 percent of the total mass of the 
galaxy. (Time in rotational periods of the initial balanced disk at r = 10 kpc.) 
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Figure 84.- Variation of the radial velocity dispersion with radius for the system in 
figure 83. (Time in rotational periods of the initial balanced! disk at r = 10 kpc.) 
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Figure 85.- Variation of the radial gravitational field with radius for the system in 
figure 83. (Time in rotational periods of the initial balanced disk at r = 10 kpc.) 
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